3. Efectos del muestreo y fenómeno de aliasing¶

3.1. El espectro de una señal discreta¶

3.1.1. ¿Qué le ocurre al espectro de una señal continua cuando la muestreamos?¶

Sin pérdida de generalidad, consideremos para este ejemplo la señal gaussiana

\[ s(t) = \exp \left ( -\frac{t^2}{2\sigma^2} \right) \]

Muestrear es equivalente a multiplicar una señal continua por un tren de impulsos, también conocido como «peineta de Dirac»

\[ \upuparrows(t) = \sum_{m=-\infty}^\infty \delta[t - m / F_s], \]

donde \(F_s\) es la frecuencia de muestreo, es decir el inverso entre la separación de los dientes de la peineta

t = np.arange(-5, 5, step=0.0001)
s = lambda t, sigma=0.5 : np.exp(-0.5*t**2/sigma**2)

Fs = Slider(start=0.2, end=5, value=1, step=.01, title="Frecuencia de muestreo")

t_dirac = np.arange(-5, 5, step=1/Fs.value)
create_figure = lambda title : Figure(plot_width=280, plot_height=280, title=title, 
                                      toolbar_location="below", x_range=(-5, 5))
p1 = create_figure('Señal Gaussiana')
p2 = create_figure('Peineta de Dirac')
p3 = create_figure('Señal muestreada')
p1.line(t, s(t), line_width=3)
for p in [p1, p2, p3]:
    p.xaxis[0].axis_label = 'Tiempo [s]'
source = ColumnDataSource(data=dict(t=t_dirac, 
                                    tops=np.ones_like(t_dirac),
                                    sd=s(t_dirac)
                                   ))

p2.vbar(x='t', top='tops', source=source, width=0.05)
p3.scatter('t', 'sd', source=source)


callback = CustomJS(args=dict(source=source, Fs=Fs), code="""
    var t_dirac = [];
    var tops = [];
    var sd = []
    for (let i = -5.0; i <= 5.0; i+=1/Fs.value) {
        t_dirac.push(i);
        tops.push(1);
        sd.push(Math.exp(-0.5*Math.pow(i, 2)/Math.pow(0.5, 2)));
    }
    source.data['t'] = t_dirac;
    source.data['tops'] = tops;
    source.data['sd'] = sd;    
    source.change.emit();
""")


Fs.js_on_change('value', callback)
show(column(Fs, row(p1, p2, p3)))

La transformada de Fourier del tren de impulsos es

\[ \mathbb{FT}[\upuparrows(t)] = F_s \sum_{m=-\infty}^\infty \delta(f - m F_s) \]

es decir otro tren de impulsos pero en frecuencia

Pueden revisar la demostración de la transformada anterior aquí

La señal muestreada es igual a \(s(t) \cdot \upuparrows(t)\)

Para obtener el espectro de la señal muestreada recordemos que multiplicar dos señales en el dominio del tiempo corresponde a convolucionar sus transformadas de Fourier en frecuencia

\[\begin{split} \begin{align} \mathbb{FT}[s(t) \cdot \upuparrows(t)] &= \mathbb{FT}[s(t)] * \mathbb{FT}[\upuparrows(t)] \nonumber \\ &= S(f) * F_s \sum_{m=-\infty}^\infty \delta(f - m F_s) \nonumber \\ &= F_s \sum_{m = -\infty}^{\infty} S(f - m F_s) \nonumber \end{align} \end{split}\]

Para entender la última igualdad revisemos la siguiente sección

3.1.2. ¿Qué significa convolucionar con un tren de impulsos?¶

Matemáticamente la convolución entre dos señales discretas (de una dimensión) se define como

\[ (f * g)[n] = \sum_{m=-\infty}^\infty f[m] g[n-m] = \sum_{m=-\infty}^\infty f[n-m] g[m] \]

Podemos imaginar que \(g\) se desplaza sobre \(f\) (o viceverza)
En cada paso el \(g\) desplazado se multiplica punto a punto con \(f\) y luego se suman
El resultado de convolucionar \(f\) con \(g\) es una nueva función

Si \(f\) es un tren de impulsos ocurrirá una «repetición» de \(g\) (o viceverza)

La siguiente animación muestra en la imagen superior un pulso cuadrado que se desplaza sobre sren de impulsos en la figura superior. La figura inferior muestra el resultado de la convolución

%%capture
fig, ax = plt.subplots(2, figsize=(7, 4), sharex=True, tight_layout=True)
t = np.arange(-4, 4, step=1e-4)

def my_signal(t, a=0, T=0.5):
    s = np.zeros_like(t)
    s[np.absolute(t-a)<T] = 1
    return s

# 2 segundos de espacio entre los dientes de la peinte
def peineta(t):
    s = np.zeros_like(t)
    s[::20000] = 1 
    return s
    
conv_s = np.convolve(my_signal(t), peineta(t), mode='same')

def update(a = 0): 
    ax[0].cla(); ax[1].cla()
    p1, p2 = my_signal(t, 0.1*a - 4), peineta(t)
    ax[0].plot(t, p1, label='señal'); 
    ax[0].plot(t, p2, label='peineta'); 
    ax[0].legend()
    ax[1].plot(t, conv_s); 
    ax[1].set_xlabel('Tiempo [s]'); 
    ax[1].scatter(0.1*a -4, np.sum(p1*p2), s=100, c='k')
    return ()
    
anim = animation.FuncAnimation(fig, update, frames=80, interval=100, blit=True)

HTML(anim.to_html5_video())

Por lo tanto: Muestrear una señal en el tiempo hace que su espectro \(S(f)\), sea cual sea, se repita infinitamente. Además el espectro se repite cada \(F_s\) [Hz]

3.1.3. ¿Cúal es el espectro de la señal Gaussiana?¶

La transformada de Fourier de

\[ s(t) = \exp \left ( -\frac{t^2}{2\sigma^2} \right) \]

es

\[ S(f) = \mathbb{FT}[s(t)] = \sqrt{2\pi}\sigma \exp \left ( -2\sigma^2 \pi^2 f^2 \right) \]

La función gaussiana es muy especial: El espectro de una gaussiana es otra gaussiana

Puedes ver la demostración de la transformada anterior aquí

En base a la gráfica siguiente notemos que

Mientras más «ancha» sea la gaussiana en el tiempo (\(\sigma\) pequeño) más «angosto» será su espectro en frecuencia
Las gaussianas son del mismo ancho cuando \(\sigma = \frac{1}{\sqrt{2\pi}}\)

x = np.arange(-5, 5, step=0.001)
sigma = Slider(start=0.1, end=2, value=1/np.sqrt(2*np.pi), step=.01, title="sigma")

source = ColumnDataSource(data=dict(x=x, 
                                    gt=np.exp(-0.5*x**2/sigma.value**2),
                                    gf=np.exp(-2*(np.pi*x*sigma.value)**2)*np.sqrt(2*np.pi)*sigma.value
                                   ))

create_figure = lambda title : Figure(plot_width=350, plot_height=280, title=title, 
                                      toolbar_location="below", x_range=(-5, 5))
p1 = create_figure('Dominio de tiempo', )
p2 = create_figure('Dominio de frecuencia')
p1.line('x', 'gt', source=source, line_width=3)
p2.line('x', 'gf', source=source, line_width=3)

callback = CustomJS(args=dict(source=source, sigma=sigma), code="""
    var x = source.data['x'];
    var gt = source.data['gt'];
    var gf = source.data['gf']
    for (var i = 0; i < x.length; i++) {
        gt[i] = Math.exp(-0.5*Math.pow(x[i]/sigma.value, 2));
        gf[i] = Math.sqrt(Math.PI*2)*Math.exp(-2*Math.pow(Math.PI*x[i]*sigma.value, 2))*sigma.value;
        
    }
    source.change.emit();
""")


sigma.js_on_change('value', callback)
show(column(sigma, row(p1, p2)))

3.1.4. ¿Cuál es el espectro de la gaussiana «discreta»?¶

Como ya vimos el espectro de una señal discreta es idéntico al de la señal continua pero «repetido» según el valor de la frecuencia de muestreo

Por ejemplo si \(\sigma=1\) y \(Fs = 2\) [Hz] tendríamos lo siguiente

def gaussiana(f, sigma):
    return np.sqrt(2*np.pi*sigma**2)*np.exp(-2*(np.pi*f*sigma)**2)

def espectro_discreto(f, sigma):
    S = np.zeros_like(f)
    for m in range(-20, 20):
        S += gaussiana(f - Fs*m, sigma)
    return S

def ventana(f, Fs):
    SW = np.zeros_like(f)
    SW[np.absolute(f) < Fs/2] = 1
    return SW

Fs = 2.
sigma = 1.
f = np.arange(-3*Fs, 3*Fs, step=1e-3)

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, gaussiana(f, sigma),  color='green', alpha=0.75,
        line_width=3, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, sigma), line_width=3, alpha=0.75,
        legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
        line_dash='dashed', legend_label='Ventana cuadrada')

p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)

Por lo tanto: Si multiplicamos el espectro discreto (azul) por una ventana cuadrada de ancho \(F_s\) (negro punteado) podemos recuperar el espectro original (verde) sin pérdidas

Esta es la base del siguiente teorema fundamental

3.2. Teorema del muestreo¶

Sea una señal continua \(s(t)\) muestreada a \(F_s\) [Hz] produciendo una señal digital \(s[n] = s(t = n/F_s)\)

Si

\[ f_{\text{max}} < \frac{F_s}{2}, \]

donde

\(f_{\text{max}}\) es la componente de frecuencia más alta de la señal
\(\frac{F_s}{2}\) es la frecuencia de Nyquist

entonces la señal analógica \(s(t)\) puede ser recuperada perfectamente a partir de sus muestras discretas \(s[n]\)

Además el valor de la señal continua reconstruida es:

\[ s(t) = \sum_{n=-\infty}^{\infty} s[n] \text{sinc}(\pi F_s (t - n /F_s) ) \]

3.2.1. ¿Qué pasa con el espectro «discreto» si no se cumple la condición anterior?¶

Asumamos que la frecuencia de muestre se mantiene \(F_s=2\) [Hz] y que \(\sigma\) disminuye

La frecuencia máxima de la gaussiana es la «última» frecuencia donde el espectro es distinto de cero

Si \(\sigma\) disminuye la \(f_{\text{max}}\) aumenta

Fs = 2.
sigma = 0.4
f = np.arange(-3*Fs, 3*Fs, step=1e-3)

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, gaussiana(f, sigma),  color='green', alpha=0.75,
        line_width=3, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, sigma), line_width=3, alpha=0.75,
        legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
        line_dash='dashed', legend_label='Ventana cuadrada')

p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)

Si asumimos que \(\sigma =1 \) se mantiene y que la frecuencia de muestreo disminuye se obtiene un efecto equivalente

Fs = 0.75
sigma = 1.
f = np.arange(-3*Fs, 3*Fs, step=1e-3)

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, gaussiana(f, sigma),  color='green', alpha=0.75,
        line_width=3, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, sigma), line_width=3, alpha=0.75,
        legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
        line_dash='dashed', legend_label='Ventana cuadrada')

p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)

Debido al solapamiento en el espectro discreto (azul) se vuelve imposible recuperar el espectro original (verde) sin alteraciones. Podríamos recuperar su parte central, pero eso significa perder información de alta frecuencia

El solapamiento espectral se llama aliasing

3.3. Aliasing¶

El espectro de una señal muestreada es periódico en \(F_s\)
Si originalmente la señal tenía componentes con frecuencias mayores a \(\frac{F_s}{2}\) se produce un «traslape» o «solapamiento» espectral
Este fenomeno se llama aliasing y los componentes traslapados se denominan aliases

Cuando existe aliasing veremos que no es posible reconstruir la señal original sin ambiguedad

3.3.1. Espectro teórico de una sinusoide¶

Consideremos por ejemplo la siguiente señal sinusoidal

\[ s(t) = \cos(2\pi f_0 t) \]

La transformada de Fourier de coseno es un impulso en \(f_0\) y otro en \(-f_0\)

\[ S(f) = \frac{1}{2} \left(\delta(f-f_0) + \delta(f+f_0) \right) \]

Puedes ver la demostración de esta transformada aquí

3.3.2. Espectro de una sinusoide discreta¶

Si muestreamos a una frecuencia \(F_s\) que sea menor a \(2 f_0\) entonces habrá traslape en el espectro discreto

Por ejemplo consideremos \(f_0 = 1.23\) [Hz] y \(F_s = 2\) [Hz]. Digamos además que observamos la señal por \(100\) [s]

¿Cómo se ve el espectro de amplitud?

import scipy.fft as sfft

f0 = 1.23
T = 100
Fs = 2 

ts = np.arange(0, T, step=1/Fs)
signal = lambda t, f : np.cos(2.0*np.pi*f*t)
SA = sfft.fftshift(np.absolute(sfft.fft(signal(ts, f0))))
freq = sfft.fftshift(sfft.fftfreq(n=len(ts), d=1/Fs))

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(freq, SA, line_width=3)
p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)

Aparece un componente en \(f = \pm 0.77\) [Hz]. Pero la señal era de \(1.23\) [Hz] ¿Puedes explicar porqué?

3.3.3. Espectro traslapado de la sinusoide discreta¶

Recordemos que el de la señal de la señal «original» se repite cada \(F_s\) [Hz]

En este caso \(F_s\) es menor que \(2 f_0\) por lo que ocurrirá traslape espectral

f = np.arange(-3*Fs, 3*Fs, step=1e-3)

def espectro_coseno(f, f0):
    S = np.zeros_like(f)
    S[np.isclose(f, -f0)] = 1
    S[np.isclose(f, f0)] = 1
    return S

def espectro_discreto(f, f0):
    S = np.zeros_like(f)
    for m in range(-20, 20):
        S += espectro_coseno(f - Fs*m, f0)
    return S

def ventana(f, Fs):
    SW = np.zeros_like(f)
    SW[np.absolute(f) < Fs/2] = 1
    return SW

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, espectro_coseno(f, f0),  color='green', alpha=0.75,
        line_width=4, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, f0), line_width=3, alpha=0.75,
        legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
        line_dash='dashed', legend_label='Ventana cuadrada')

p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)

El cuadrado negro con linea punteada representa el rango

\[ \left[-\frac{F_s}{2}, \frac{F_s}{2} \right] \]

es decir lo que «podemos ver» cuando usamos la transformada de Fourier (es equivalente a la figura anterior)

3.3.4. «Aliases» de la señal¶

Como vimos en este caso la frecuencia de la señal está fuera del rango que podemos estudiar

Además por culpa del solapamiento espectral aparecerán «aliases»

La frecuencia de los aliases sigue la fórmula

\(f_a = f_0 + m F_s\) [Hz] donde \(m\) es un número natural
\(f_a = m F_s - f_0\) [Hz] donde \(m\) es un número natural

Uno de estos aliases es \(f_a = Fs - F0 = 0.773\) [Hz], la frecuencia que vimos originalmente

T = 5
t = np.arange(0, T, step=1e-4)
ts = np.arange(0, T, step=1/Fs)
signal = lambda t, f : np.cos(2.0*np.pi*f*t)

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")

p1.line(t, signal(t, f0),  color='green', alpha=0.75, 
        line_width=3, legend_label='Señal original')
p1.scatter(ts, signal(ts, f0), size=10, legend_label='Señal muestreada')
p1.line(t, signal(t, Fs - f0), line_width=3, alpha=0.75, color='black',
        line_dash='dashed', legend_label='Alias de la señal')

p1.xaxis[0].axis_label = 'Tiempo [s]'
show(p1)

El «alias» se hace pasar por la señal original. No podemos distinguir cual de las dos se uso para generar la señal muestreada (puntos azules). Ya no podemos reconstruir sin ambiguedad la señal

3.3.5. ¿Cómo eliminamos el aliasing?¶

Necesitamos que todas las frecuencias del espectro sean menores que \(\frac{Fs}{2}\)

Podemos lograr esto mediante

Filtrado: Eliminar las frecuencias mayores a \(\frac{F_s}{2}\) (Próxima unidad)
Aumentar \(F_s\) tal que sea dos veces mayor que la frecuencia máxima de interés, no siempre es fácil saber cuál será su valor a priori

3.4. Apéndice: Principio o Teorema de incertidumbre¶

El principio de incertidumbre de Heisenberg nos dice que la precisión (certeza) con que medimos la posición de una particula es inversamente proporcional a la precisión con que medimos su momentum lineal:

\[ \Delta x \Delta p \geq \frac{h}{4\pi}, \]

donde \(h\) es la constante de Planck

En señales existe un principio análogo: No podemos especificar con infinita precisión la localización temporal y frecuencial de una señal al mismo tiempo.

Denis Gabor (1946) fue el primero en darse cuenta de que el principio de incertidumbre aplica para señales.

Su teorema dice que para una señal con energía finita

\[ E = \int |s(t)|^2 dt \]

con valor medio temporal

\[ \langle t \rangle = \frac{1}{E} \int t |s(t)|^2 dt, \]

y varianza temporal

\[ (\Delta t)^2 = \frac{1}{E} \int (t - \langle t \rangle)^2 |s(t)|^2 dt, \]

cuya transformada de Fourier \(\mathbb{FT}[s(t)] = S(\omega)\) tiene un valor medio en frecuencia

\[ \langle \omega \rangle = \frac{1}{E} \int (\omega - \langle \omega \rangle) |S (\omega)|^2 d \omega \]

y varianza frecuencial

\[ (\Delta \omega)^2 = \frac{1}{E} \int (\omega - \langle \omega \rangle)^2 |S(\omega)|^2 d\omega \]

se cumple que

\[ \Delta t \Delta \omega \geq \frac{1}{2}, \]

es decir \(\Delta t\) y \(\Delta \omega\) no pueden ser arbitrariamente pequeños

Tal como vimos antes el ancho temporal y el ancho frecuencial están inversamente correlacionados

INFO183 Análisis de sistemas lineales